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Transformed snapshot interpolation 


G. Welper* 


Abstract 


Functions with jumps and kinks typically arising from parameter de¬ 
pendent or stochastic hyperbolic PDEs are notoriously difficult to approx¬ 
imate. If the jump location in physical space is parameter dependent or 
random, standard approximation techniques like reduced basis methods, 
PODs, polynomial chaos, etc. are known to yield poor convergence rates. 

In order to improve these rates, we propose a new approximation scheme. 

As reduced basis methods, it relies on snapshots for the reconstruction 
of parameter dependent functions so that it is efficiently applicable in a 
PDE context. However, we allow a transformation of the physical coordi¬ 
nates before the use of a snapshot in the reconstruction, which allows to 
realign the moving discontinuities and yields high convergence rates. The 
transforms are automatically computed by minimizing a training error. 

In order to show feasibility of this approach it is tested by Id and 2d 
numerical experiments. 

Keywords: Parametric PDEs, reduced order modelling, shocks, transforma¬ 
tions, interpolation, convergence rates, stability 
AMS subject classifications: 41A46, 41A25, 35L67, 65M15 

1 Introduction 

A cornerstone of reduced order modelling, stochastic PDEs and uncertainty 
quantification, is the efficient approximation of high dimensional PDE solu¬ 
tions u(x, fi) depending on physical variables x £ Q and parametric or random 
variables [i £ V C M. d . Many contemporary approximation techniques like e.g. 
reduced basis methods [2H1221E3] , POD [271 US QH], Karhunen Loeve expansion 
m or polynomial chaos [2HI ESS G2] build upon a reconstruction by a truncated 
sum 



( 1 ) 


where the choice and computation of Cj(/z) and ipi{ x ) depends on the specific 
method at hand: For reduced basis methods ipi(x ) = u(x,Hi) are snapshots 
and the Ci(/z) are computed by a Galerkin projection. For POD and Karhunen 
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Loeve expansions one minimizes the error between u(x,p) and any truncated 
representation of the form |l]) and in case of polynomial chaos the functions Ci(p) 
are orthogonal polynomials. Borrowing from the tensor community, we refer to 
0 as a polyadic decomposition and denote methods based on it by polyadic 
decomposition based methods. The success of all these methods relies on the 
fact that for many problems one can truncate this sum to a few summands only 
for the price of a very small error. 

However, this regularity assumption is not always true. An important class of 
problems are functions u(x , p) that have parameter dependent or random jumps 
or kinks arising e.g. in parametric or stochastic hyperbolic PDEs. Polyadic 
decomposition based methods are expected to perform poorly for these type 
of problems. In fact in Appendix [A] we consider a counterexample for which 
no polyadic decomposition based method can achieve a convergence rate higher 
than one with respect to the number of summands in a polyadic decomposition. 
See also [T2] for a survey in case of uncertainty quantification. 

There are relatively few methods in the literature ®m\m that directly 
address this poor performance for parameter dependent jumps and kinks. In¬ 
stead, much of the work does use polyadic decompositions and focuses on dif¬ 
ferent problems arising in the context of reduced order modelling of paramet¬ 
ric hyperbolic PDEs and singularly perturbed problems: Solving the PDE di¬ 
rectly in a reduced basis, online/offline decompositions and error estimators, see 

nn udied me. 

The goal of this paper is the construction of an alternative approximation 
method to replace standard polyadic decompositions in order to achieve higher 
convergence rates for functions u(x,p) with parameter dependent jumps and 
kinks. In addition, the method relies on snapshots and optionally error esti¬ 
mators as input data so that it can be used efficiently and non-intrusively with 
existing PDE solvers. Somewhat similar to £22], we allow a transformation 
<p(p,rj) : Q — > D of the physical domain before we use a snapshot u(x,r)) in the 
reconstruction of u(x,p), i.e. 

u{x,p)^ ^ 
riev„ 

where r] is in some finite set of parameters V n . The purpose of the additional 
transform is an alignment of the discontinuities of u(x,i 7 ) with the ones of 
u{x,p). As a result the discontinuities are “invisible” in parameter direction 
so that very few summands yield accurate approximations. More rigorously, 
we prove a high order error estimate that does not depend on the regularity of 
u{x , p) itself, but on the regularity of the modified snapshots after alignment 
which is considerably higher for many practical problems. In addition, because 
exact alignment is rarely possible in practice, we also take perturbation results 
into account. Similar to greedy methods for the construction of reduced bases, or 
neural networks, the transform </> is computed by minimizing the approximation 
error on a training sample of snapshots. Although this might seem prohibitively 
complicated, in Section[5]we discuss some preliminary arguments to avoid being 
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trapped in local minima and in Section [g] some 2d numerical experiments are 
provided where simple subgradient methods provide good results. 

The outlined approximation scheme allows for various realizations with re¬ 
gard to the choices of the coefficients or the inner transforms 

Because the main objective of the paper is a proof of principle that one can 
approximate functions with parameter dependent jumps and kinks with high 
order from snapshots alone, we usually vote for the simplest possible choices 
and leave more sophisticated variants for future research. 

The paper is organized as follows: In Section [2] we present the main approxi¬ 
mation scheme and a basic error estimate. Then, in Section[3]we prove stability 
results and consider approximations of the inner transform. Afterwards we 
turn to the actual construction of </>(/x, 77). We first discuss some characteristic 
based approaches and their drawbacks in Section [4] and then the optimization 
of </>(/x, 77 ) by training errors in Section [ 5 J Finally Section [6] provides some Id 
and 2d numerical experiments. For the sake of completeness, in Appendix |A| we 
consider a counterexample for which no polyadic decomposition based method 
can achieve high order convergence rates. 


2 Transformed snapshot interpolation 


In order to motivate the new approximation scheme, we consider the following 
prototype example throughout this section 


u(x , /i) := ip 


0.4 +11 




-1 < x < - \ 

else 


( 2 ) 


where ip is the standard modifier cut off at x = —1/2. This is not necessarily a 
solution of a PDE, but has the main features we are interested in: a discontinuity 
that is moving with the parameter. An example for a polyadic decomposition 
based approximation can be seen in Figure |T] where we recover u(x , /x) from 
three snapshots by a polynomial interpolation 




where V n C V are some interpolation points and £ v , rj £ V n are the corre¬ 
sponding Lagrange interpolation polynomials. We see the typical “staircasing 
behaviour” which significantly deteriorates the solutions quality. Although our 
choice of the polyadic decomposition is perhaps overly simplistic, reduced basis 
methods and other more sophisticated schemes suffer from the same problem. 

Unlike this superposition of snapshots resulting in the “staircasing” phenom¬ 
ena, it seems much more intuitive to compute one snapshot u(x 1 rj) and recover 
the function u{x , /x) for a different /x by stretching this snapshot such that the 
left end of the support is fixed and the jump locations match. To state this 
intuition in mathematical terms, “stretching” one function u(x,ij) to match a 
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Figure 1: Left: Snapshots of the parametric function © for parameters fx = 
0.5,0.85,1.2. Right: Exact solution (red) and polynomial interpolation (blue, 
dashed) for \x = 1.0. 


second one u(x,fx) essentially boils down to a transform </>(/x, 77) : £2 — > fl of the 
physical variables so that we obtain the approximation 

u(x,/x ) « u((j>(fx,r])(x),r]). 

In general, even for optimal choices of (j), we cannot obtain arbitrarily good 
approximation errors in this way. To obtain convergence, we therefore require 
in addition that 77 is close to /i which yields the following simple approximation 
scheme: First, we choose a finite subset V n C V of the parameter domain V 
and compute the snapshots u(x,rj) for r) £ V n . Then, given a new /i € P, we 
find the ? 7 M £ V n closest to fx and approximate 

u(x , fx) « u((/>(n, r]n){x), Vn)- (3) 

Besides the snapshots themselves this requires the knowledge of the transforms 
(x, fx) —> <p(fx , 77) ( x) for finitely many 77 £ V n - Thus instead of approximating one 
single function depending on x and /i we have to find many of them! However, 
whereas polyadic decompositions perform poorly for u(x,n ), they often yield 
good results for the transforms <(>(p., rf) : Their smoothness with respect to /i 
depends on the smoothness of the jump or kink location with respect to the 
parameter and not the smoothness of u{x , /z) itself. For example ([2]) the jump 
location is j (/x) = | + ^/x so that a linear transform 

= x - j(/x) + j(rj) 

is sufficient to align the jumps. As shown in Figure [2j this transform does not 
align the left end of the supports, but because r\ is close to /z this is good enough 
as we see below. For more complicated problems the transform <(>(/z, rf) is not 
explicitly known and we have to find efficient ways to compute it from the given 
data. We postpone this issue to Section [5] and assume for the remainder of this 
section that 4>{ix,rf)(x) is given to us. 
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Figure 2: Left: Transformed snapshots of the parametric function (|2j) for param¬ 
eters /i = 1.0 and r\ = 0.5, 0.85,1.2. Right: Exact solution (red) and transformed 
snapshot interpolation (blue, dashed) for 77 = 1 . 0 . 


To assess the approximation error, we observe that our scheme is a 
piecewise constant approximation of the transformed snapshots 

v^{x,rj) := u(<t>(n,rj)(x),ri) 

with respect to 77 at the point 77 = 77. Thus, we obtain the error estimate 

IMz, AO ~ v„(x, Vu.)\\l p = Oin- 1 ) (4) 

where n is the number of snapshots, provided that v^x, rf) is differentiable with 
respect to 77. This is achieved by the inner transform </>(//, 77): Whereas the 
original snapshots 77(2 , 77) have jumps in parameter direction, the transformed 
snapshots have jumps in fixed locations independent of 77 resulting in a smooth 
dependence of 7^(2, 77) on 77. Because 7/7(77,77) is supposed to align the disconti¬ 
nuities and kinks of u(-, 77) and u(-, 77), it is natural to require that 

ix){x) = x (5) 

which yields 77(2,77) = v^x, 77). With Q, it follows that 

\\u(x,n) - v ll (x,r] ll )\\ Lp = C>(n _1 ) 

so that our approximation scheme achieves first order convergence. 

In comparison, due to lacking smoothness for standard piecewise constant 
approximations 

77(2,77) « u(x, 77^) (6) 

we expect convergence rates of h 1 ' p for spacial errors in L p . Thus, depending 
on the norm, the inner transform yields a gain in the convergence order for 
p < 1 or none at all for p = 1. However, the major impediment is no longer 
a lack of regularity but the low order convergence of the piecewise constant 
approximation of 77 — > 77^(2,77). Therefore, we replace it by a higher order 
scheme. For simplicity, we confine ourselves to a simple polynomial interpolation 
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and leave more sophisticated choices for future research. Thus for interpolation 
points V n C V and corresponding Lagrange basis polynomials £ v we define the 
transformed snapshot interpolation by 

u(x,fj,) « u n (x,ii) := ^2 ( 7 ) 

r)EV n 


The input data for this reconstruction is identical to the previous piecewise 
constant case: We need \V n \ snapshots and \V n \ transforms (x,p) —> (j>(p,r])(x), 
r/ € V n . Only the reconstruction formula has been changed to a higher order 
interpolation. In order to state an error estimate, let P" be the span of the 
Lagrange basis polynomials and recall that the Lebesgue constant is the norm 
of the polynomial interpolation operator in the sup-norm which is given by 

n 

A„ := sup V \£ v (p)\. (8) 


We obtain the following error estimate. 


Proposition 2 . 1 . Assume u n (x,p) is defined by the transformed snapshot in¬ 
terpolation 0 . Then for all p £ V the error is bounded by 


I \u(‘,ll) - U n (-,n)\\ Ll < A r 


mf n K, ( 0-p| ^ 


Proof. The proof follows directly from u(x,p) = v^{x,p) and standard interpo¬ 
lation estimates applied to rj —> v^x^rf). □ 

For this proposition as well as the remainder of this paper we choose the 
Li-norm to measure errors because it is the most common choice for hyperbolic 
PDEs. Also note that the given result is just one option of the various estimates 
for polynomial interpolation. For example, if we assume analytic dependence 
of v fl (x,r]) on ?y and use Chebyshev nodes for a one dimensional parameter, we 
can achieve exponential convergence rates. The most important observation, 
however, is that the estimate does not involve any regularity assumption of 
u(x,p) itself. Instead it relies on the regularity of Vn(x,r]) with respect to 77 
which can be considerably better. 

The results for example 0 are shown in the right picture in Figure [2j We 
see a very accurate approximation of the jump, however the approximation 
quality around the left end of the support of u(x, /i) is slightly worse than for 
the original interpolation in Figure [Tj The reason is that the left end of the 
support is parameter dependent after the transform so that v^x, p) is no longer 
analytic in 77, however infinitely differentiable. Therefore the loss we suffer at 
this point is of orders of magnitude less than the staircasing behaviour around 
the jump of the simple interpolation in Figure [l] 

In summary, instead of approximating the non-smooth function u(x, /i) di¬ 
rectly, for every target p £ V we construct a new smooth function (x, rf) — > 
Vfi{x,rj) and approximate this function instead. The interpolation condition 0 
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guarantees that u{x,fi) = v lx (x, rf) so that this yields accurate approximations 
of 77(2,77) itself as depicted in Figure [3j In addition, for our preliminary simple 
linear interpolation of 77^(2,77) this allows an offline/online decomposition: in 
an offline phase, we compute the snapshots u(x , 77) as well as the transforms 
ct>([i,ri)(x ) at the interpolation nodes 77 (see Section [ 5 ] below). Then in an online 
phase we can efficiently approximate 77(2, /z) for any /z £ V by the transformed 
snapshot interpolation (J7|. 



Figure 3: The three dimensional coordinate axes indicate the linear space of 
x-dependent functions so that each point of the red line represents a snapshot 
2 — > u( 2,77) for some parameter 77. The kinks indicate that in our case this 
solution manifold is not smooth with respect to 77, however note that in reality 
it is non-smooth for every parameter. The blue dashed line indicates the more 
smooth transformed snapshots 7^(2,77). 


3 Stability 


Of course high order smoothness of 7^(2, rf) with respect to 77 needed for high ap¬ 


proximation orders in Proposition 2.1 requires that jumps and kinks are exactly 
aligned. However, for any finite approximation of the inner transform </>, this is 
rarely possible. Therefore, we next consider two perturbation results, that allow 
us to bound the error while taking approximation errors of the inner transform 
into account. The first one, Lemma [3T] relies on a measure theoretic argument 
and allows rather general transforms including ones with kinks as found in e.g. 
finite element discretizations. The second one, Corollary 3.2 avoids measure 


theory, but requires the inner transforms <)>(/z, rf) to be diffeomorphisms. 

In the following, let tp(fj,,i])(x) be a perturbation of <t!)(/z, 77)(2) . To simplify 
the arguments below, for the time being, we forget about the parameter depen¬ 
dence and consider two transforms cf> : fl —> fl instead. If we assume that 
each point <f>(x) can be connected to the point 73(2) along a curve <h s (2) for s in 
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the interval [0,1], we can rewrite the perturbation by the fundamental theorem 
of line integrals 

(u o cf>) (x) — {u o ip)(x) = f u'(<b s (x))d s $ s (x)dt, 

Jo 

so that it remains to estimate the right hand side. The map (x, s) —> $ s (x) can 
be regarded as a function from f2 x [0,1] —> f2 such that 

<b 0 (x) = <j>{ x) $ 1 (x) = <p(x). (9) 

which is a homotopy between (j> and ip if it is continuous in addition. Further¬ 
more, let A be the Lebesgue measure and A the Lebesgue u-algebra on 17 and 
let TJA denote the pushforward measure defined by 

$:A(A) = A(($ S )- 1 (A)) for all A £ A. 

Then we have the following lemma. 

Lemma 3.1. Assume that u £ BV(fl) and <I> S , 0 < s < 1 given by ([9]) is 
measurable and differentiable with respect to s such that 

$*A(Al) < cX(A) for all A £ A and 0 < s < 1 (10) 

and 

sup |9 s 4> s (a;)| < C\\<j) - <p\\ Loo (n) (H) 

0<s<l 

for constants c, C > 0. Then we have 

\\uo<j>-uoip\\ Ll ^ < cC'||u|| w( n)||^-¥’|U 00 (n)- (12) 

Let us discuss the main assumptions before we prove the proposition. If the 
speed |9 s <E> s (:r)| of each curve s —> $ s (ir) is quasi uniform, i.e. equivalent to a 
constant S for all x and s, we have 

||9 s $ s ||L 00 (nx[o.i]) ~ [ |<9 s $ s (x)|d s=:l(x) 

Jo 

where l(x) is the length of the curve connecting <f(x) to ip(x). In that case 
assumption © states that, up to a constant, the length of each curve <b s (x) is 
bounded by the distance \<j)(x) — <p{x)\ of its endpoints. 

In case the domain 17 is convex, a simple choice of the curves d> s (x) are the 
convex combinations of the end points: 


>F s (x) = (1 — s)(j>(x) + sip(x). 


In that case, we have <9 s <l> s (x) = <p(x) — <j>(x) so that condition (11) is satisfied. 

In order to justify the second assumption (10), let us consider the following 
scenario: Assume that <j>{x) = xq £ 17 and (p(x) = x± £ O map all of 17 to 




single points. Furthermore let u be a piecewise constant function with a jump 
so that xo and x\ are on different sides of this jump. On the one hand we obtain 
||«o <fi — u o = ||li(rco) ^ u(xi)||i 1 (f 2 ) = where h is the hight of the 

jump. On the other hand we have ||0 — ¥’||L 00 (n) = l^o —Xi| which can be made 
arbitrary small by suitable choices of Xg and X\ on each side of the jump. Thus 
the main statement of the proposition is violated. This counterexample 
relies on the fact that both transforms concentrate all weight in a single point 


such that A({xi}) = A(fi), i = 0,1 which is ruled out by assumption (101. 

Finally, we assume that the outer function u £ BV(fi) is of bounded vari¬ 
ation. This allows jumps and kinks and is one of the most common norms for 
stability results of hyperbolic PDEs. 


Proof of Lemma \3.1\ For the time being, let us assume that u £ C 1 (fl). Ap¬ 
plying the fundamental theorem for line integrals, we obtain 

(it o — (u o ip)(x) = ( ix / (3» s (£c))«9 s <P> s (a^) dt 

Jo 


Thus, we have 


\\uo(/)-uoip\\ Llin) = f ( it'($ s (a;))(9 s $ s (a;)ds 

Jn Jo 

< ( f |u , ( < h s (a;))||3 s $ s (a;)| dsdx 

Jn Jo 

< sup |9 s $ s (a;)| f I |it'($ s (a;))| dsdx 

o<s<i Jn Jo 

< C\\(/) - ^||i oo( n ) [ f |u / ($ s (a;))|dxds 

Jo Jn 

Using the pushforward measure <bJA” +1 and its bound (101 we conclude that 

f |u'($ s (a:))|dx = [ \u'(y)\ d$*A(:r) < c f |u'(y)|dx (13) 

Jo, Jq Jq 

Combining the last two estimates and using that /),' ds = 1 yields 

\\uo<j)-uoip\\ Ll ( n) < cCj|<^>-<p|| Loo(n) / \u'(y)\dy, 

Jn 


which is equivalent to the estimate (12) we wish to prove. 

Finally, we extend the estimate to all u £ BV( fl) by using a density argu¬ 
ment. To this end note that for all e > 0 there is a u e £ C 1 (C2) such that 


\u-Ue\\ Ll (n) < e 


l^eIILi(o) < IM|sv(fi) +e 
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Thus, to apply a density argument, is suffices to bound \\uo(f> — u e o^)\\ Ll /Q^ and 
||u o tp — u e o Analogously to (13) we obtain 


\\uo(j)-u e o(j)\\ Liin) = / |u($°(a;))-u e ($ 0 (a;))|dx 
Jn 

= [ I u(y) -u e (y)|d$°A(j/) 

Jn 

<c \u{y) - u e (y)\dy 
J n 

< c\\u - U e I|x, 1 (f2) 

The bound for ||u o <p — u e o follows analogously which completes the 

proof. □ 

If the transforms $ s (x) can be chosen to be diffeomorphisms, the pushfor- 
ward measure is explicitly given by the usual transformation rule 


*t\(A)= [ |detD a: ($ s )- 1 (x)|dA(x) (14) 

J A 

so that we obtain the following corollary. 

Corollary 3.2. Assume that u £ BV(fl) and that <L S , 0 < s < 1 given by § 
are diffeomorphisms for fixed s and differentiable with respect to s such that 

| det -D x (<I> s ) _ 1 (a:)| < c for all A £ A and 0 < s < 1 

and 

sup |<9 s $ s (x)| < C\\(j> - <p|| Loo (n) 

0<s<l 

xEfl 

for constants c,C > 0. Then we have 


l|w ° </> - u o ^|| Ll (n) < cC\\u\\ B v(n)\\(l> - vh^tn)- 
Proof. We just have to show the bounds ITol) of the pushforward. By its explicit 


formula (141 we have 

$*A(A) = [ | det-D x ($ s ) _ 1 (a;)| dA(x) < cA(A) 

J A 


for all A £ A and 0 < s < 1 so that the corollary follows from Lemma 3.1 


□ 


Let us now consider the transformed snapshot interpolation 0 again. As¬ 
sume that there is a transform 77 ) (x) that aligns the jumps and kinks exactly 
so that we obtain high convergence rates in Proposition |2.1| In general, we have 
to find a finite approximation to this exact transform, say <p m (y,,r])(x). Note 
that according to ([7| we only need to know this function for the \P n \ nodes 
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rj £ V n , so that we have to approximate V n functions depending on x £ Q 
and a parameter p £ V. Of course this is exactly the same problem as approxi¬ 
mating a function u(x, p) which is our initial problem, however, the regularity 
of 4 > can be much more favorable as we have seen in the introduction 
in Section [2] or as we will see in Section [4] below. Therefore, we can apply a 
more classical polyadic decomposition based approach to find an approxima¬ 
tion <f> m (p 1 r])(x) of the inner transform. Replacing the exact transform by the 
approximate one in the transformed snapshot interpolation yields 

u(x, p) « Un,m(x,n) ~ ^ t- n (p)u(<f> m (/i, T]) (x) , T}). (15) 

Combining the error estimate of Proposition |2.1| with the perturbation result 
Lemma [XT] we arrive at the following Proposition. 

Proposition 3.3. Assume thatu £ BV{fi.) and that there are curves $ s (p, r])(x), 
0 < s < 1 for x £ fl, p £ V, p £ V n measurable and differentiable with respect 
to s such that 


and 

and 


$(M, V)°(x) = </>(p, rf){x) ${p, v) 1 ( x ) = <M/b rf){x). 

$(//, J?)*A(A) < c\(A) for all A £ A and 0 < s < 1 


sup \d s ^{p,r/) s (x)\ <C\\(t>(fi,r])-(l> m {n,r])\\ Laom 

0<s<l 


for constants c, C > 0. Furthermore let u n ^ m (x, p) be defined by the transformed 
snapshot interpolation (15). Then for all p £V we have the error estimate 

\\u(-,p) - u nt m(-,p)\\ Ll < A„ inf \v(-,p)-p\ 

per " 

+ cCA„ max ||u(-, v)\\bv(Q) max || <f>{p,rj) - <t> m {p, 

r)€Vn 

where A„ is the Lebesgue constant ([8]). 

Proof. We have 

||u(-, / u)-u„ )TO (-,/i)|| Ll ( Q ) < \\u(-, p)-u n (-, p)\\ Ll ( n) + \\u n (-, p)-u njm (-, p)\\ Ll ( Q ) 

With Proposition |2.1| the first term can be estimated by 


\\u(-,p)-u n (-,p)\\ Ll <A n inf \v(-,p)-p | 

per n 

In order to estimate the second term, using the definition of the Lebesgue con¬ 
stant and Lemma 13.1 1 we obtain 


\u n (-,p) - «7i,m(•,/*)Hz,!(fi) < X! u (<t>(l J '’V)( x ),V) ~ u(<j>m(p,rj){x),Tj) 

V&Vn 


L i(fi) 


< A„ max \\u((j)(p, r])(x), rj) - u(<j) m (p,v)( x ),v)\\L 1 (Q) 

r}^V n 

< cCA n max ||w(-, ??)||w(fi) max \\</>(p,r]) - ^llwn) 

rjEVn r)€-Vn 
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Combining all three estimates completes the proof. 


0 


If the p, dependence of rj) is smooth, we can use a polyadic decomposition 
for its approximation. Although there are much more sophisticated methods, 
possibly the simplest choice is a linear interpolation 

<l>m(n,v)(x) = ( 16 ) 

v&P m 


where t v are Lagrange basis polynomials with respect to nodes in so me f inite set 


Vm- With this inner approximation, the error bound of Proposition 3.3 depends 


of the smoothness of the transformed snapshot v^x, rj) with respect to r) and of 
the transforms (x,n) —> (j)(fi,r])(x) with respect to ji. If both dependencies are 
analytic, for suitable interpolation points the error decays exponentially. 


4 Inner transforms by characteristics 

We still have to choose an inner transform </>(/x, rj) such that the transformed 
snapshots v ft (x. rj) are as smooth in ij as possible. One obvious idea that comes 
to mind is to somehow make use of characteristics. In this section, we discuss 
some problems that arise from that approach for the Riemann problem for 
Burgers’ equation. In our example, the parameter is the hight of the jump in 
the initial condition which yields the parametric PDE 



u(x, 0) = g^x) = | o x > 0 for 1 = °' 


In addition, we assume that fi > 0 so that the solution has a shock along the 
curve i/xi. To write down an explicit solution formula, let y(x,t) be the origin 
(at time t = 0) of the characteristic passing through the point (x, t). It is easily 
seen to be 


X/xOM) 


x — /A 
x 


X < hut 
f 

X > 2 /A■ 


(17) 


Because u(x, t, fi) is constant along characteristics, we obtain 


u(x,t,fj ,) = gp{Xn(x,t))- 


(18) 


The previous discussion aside, a simple idea for an approximation scheme is to 
encode or approximate g^{x) and the characteristics \^{x,t) and then use the 
exact solution formula (18) to reconstruct u(x,t,n). In spirit this is similar to 
our original idea ([3]) where the snapshots u(-,rf) are replaced by g v {-) and the 
transform 0(/x, 77 ) by the characteristic However, from the explicit formula 
(17) for y /t , we see that this function has a parameter dependent jump. Thus, 
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in general, we have to face the same difficulties for approximating the param¬ 
eter dependent characteristic (a i,t,y) as for the original solution 

(x, t , y) —> u(x, t, y) so that there is no progress with respect to this issue. 

If we want to use characteristics to define the inner transform <f>(y, g) of the 
transformed snapshot interpolation, the problems are even more complicated. 


As for (18) we can follow the characteristics backward in time, but because we 


to not evaluate the initial condition g M but a snapshot u(x,t : g), we then follow 
the characteristics forward in time with a different parameter. To this end, let 
ifin (y, t) be the position of the characteristic at time t, starting at the initial 
position y at time t = 0: 





-\gt<y < \yt 
\l-it < y. 


(19) 


Then we can transform one solution for parameter g into a solution for parameter 
9 by 

u{x, t, n) = (20) 

g g 

However, this formula is only correct for y>g. The reason is that the interval 
of points at t = 0 that eventually end up in the shock at time t is strictly larger 
for larger parameters. Thus, for y < g there is a interval I around the shock 
location of y for which x /i (I, t) is mapped into the shock location \g by the 
forward characteristic Thus, the right hand side of (20) has only one single 
value in the interval I or is undefined whereas the left hand side has to different 
values and the formula is thus not correct. 

Nonetheless, for /i > g, we can define the transform 


<K/b y)(x, t) = <P v (Xii{x, t), t ) 
so that by (p0|) the transformed snapshot is 


( 21 ) 


v„.(x,t,g) = u(<j>(n,g)(x,t),g) = u{p v {Xti{x,t),t),t,g) = -u(x,t,n) 

I - 1 

which is clearly smooth and in fact even linear in g. Using y > g it is easy to 
verify that g) is 


<j>(y,g)(x,t) = 


x — {y — g)t x < \yt 
x x > 


Note that for our approximation scheme (J7| we need to know the function 
(x,y) —> <p(y,g)(x) for finitely many g £ V. Again, this function has a y 
dependent jump so that its approximation poses the same difficulties already 
encountered for u(x,t,y) itself. 


However, we are not obliged to use the transform (21) based on characteris¬ 
tics. By noting that the shock location is \yt, simply shifting the whole solution 
in x-direction by 

(/>(y,g)(x,t) = x + \(y-g)t (22) 
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aligns the shocks, i.e. the transformed snapshot w(</>(/i, 77 )(ar, t),t, if) has its shock 
in the location which is the shock location for parameter yt. In addition the 
interpolation condition ([5| is obviously satisfied. Note that by ( |l7| ) and (181 
the parametric solution is 

u(x,t,v.) = {l 


so that the transformed snapshot becomes 


v»(x 1 t,r 1 ) = \ l = \(x,t,») 


which is the same as for the characteristics based transform, but now for all 
£ V . Recall that the error estimate of Proposition |2.1| just requires 
smoothness with respect to r) which is obviously the case. However, opposed to 
the characteristic construction also <p(fj,,rj)(x) is smooth in p, so that polyadic 
decomposition based methods yield accurate approximations of <f> at low cost. 


5 Optimizing the interpolation error 

5.1 Generalized gradients 

We still need to find a realistic way to actually compute the inner transform 
<j>(li, if). Similar to the construction of reduced bases, PODs or neural networks, 
we aim at finding an inner transform <f> that minimizes the approximation error. 
To this end, we measure the error in the sup-norm with respect to the parameter 
which is typical for reduced basis methods but not mandatory. It follows that 
the overall error is given by 


:= sup (</>), 

IJ.&V 


(23) 


where 

cr M (</>) := || u{-,n) - u n {-,n](t >)|| Ll(n) 

is the error for one fixed parameter. To make the dependence on the inner 
transform more explicit, in this section we denote the transformed snapshot 
interpolation |7| by u n (x, /x; (j>) = u n (x,/j,). In practice it is not possible to 
minimize the error a-p (c f> ) directly because it would require the knowledge of all 
functions it(-, p) for all parameters /j. £ V. To this end, we only assume to know 
the errors a of a finite training sample /i £ T C V so that the overall error 
(23) is replaced by the training error 


o- t{4>) '■= sup O>(0). 

AiGT 


(24) 


Although surrogates for the training error are available for some singularly per¬ 
turbed problems [ZD Eol 0 E] we omit these in favor of future research. Instead, 
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we resort to an explicit knowledge of some training snapshots p, £ T 

in addition to the snapshots that are used for the reconstruction ([7]) itself. In 
contrast to the reduced basis method this severely limits the size of the training 
sample. Nevertheless, in Section[6]we consider examples which yield good results 
with roughly twice as many training snapshots than reconstruction snapshots, 
so that the additional burden of the training samples is reasonable. 

Because we are explicitly interested in non-smooth functions u, the error 
crj -(</>) is a non-trivial objective function to minimize. The next proposition 
shows that despite possible jumps of u the error is Lipschitz continuous, nonethe¬ 
less. The assumptions of this proposition are essentially the same as for Lemma 
|3.1| ancl are commented right after it. 

Proposition 5.1. Assume thatu £ BV(Q) and that there are curves <b s (p,, r])(x), 
0 < s < 1 for x £ fl, fi £ V, r/ £ V n measurable and differentiable with respect 
to s such that 

®°(p,v)(x) = = ¥>(/L»7)(aO- (25) 

and 

r/)*\(A) < cA(A) for all A £ A and 0 < s < 1 

and 

sup \d s ${/i,r]y(x)\ < C\\<j>(p, rj) - 7)|U 00 (fi) 

0<s<l 

for constants c,C > 0. Then we have 

\a r {<f) - a r (<p)\ < cCA n sup ||u(-, v)\\bv(Q) sup ||</>(/r, rf - ip(p, v)\\ Loo <n) ■ 

rieVn n&v 

U&'Pr, 

(26) 


Proof. Note that the triangle inequality implies that 


sup a^{4>) - sup a^p) 
/^GT - /^gt - 


< sup |0-„(0) -<r M (v>)| (27) 

ue T 


Wr(<t>) - °t{<p) | = 

so that it is sufficient to bound |cr At (^>) — a^fp)\. To this end note that 
WM) - IUi(n) ^ II u(-,n) - u n {-,n\ip) || Ll( n) 

L i(ii) 

H<K/b V)(x), rf) - u(ip(iJ,, v)(x),v)\ 


< 


< 


T/eVn 


L i(n) 


< -u{<p{/J,,r))(x),7i)\\ Li(n) 

riGV „ 
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Due to the given assumptions, we can now apply Lemma 3.1 to conclude that 


WM) -c^O)] < cC - v{v,v) || Loo( n) 

v&v n 

<cCA n sup |K-,??)||w(n) sup \\<f>(n, rj) - <p(fi,v)\\ Lao (n) 

T|6P„ 


With (27) this yields claimed estimate (26). 


□ 


In order to optimize the training error crj- (</>), we search for a minimizer in 
a set of candidate transforms (f) £ $ C X in a Banach space X. The space of 
continuous functions C(V xPxll) with the additional restrictions from Propo¬ 


sition 5.1 seems to be a reasonable choice for $ because in the last proposition 
the transform error is measured in the supremum norm. In general this objec¬ 
tive function is not differentiable so that we cannot rely on standard gradient 
based optimizers. However, because cttW is Lipschitz continuous according to 
the last proposition, we can use optimization methods from non-smooth opti¬ 
mization [Hill relying on the generalized Clarke gradient [3J. To this end, for 
a direction v £ X 7 we first define the generalized directional derivative 


’(</>; v) := limsup 
</?—>•</>; h.4-0 


+ hv) — cr(<p) 


where we suppress the additional T subscript of cr for simplicity. Note that 
this limit is well defined because cr is Lipschitz continuous. In order to define a 
gradient from these directional derivatives, recall that in the differentiable case 
one can define the gradient Vrr £ X* variationally by 


d v a = (Vcr, v ), for all v £ X 

where X* is the dual space of A' and (•,•) the corresponding dual pairing. 
Likewise, in the Lipschitz continuous case we define the generalized gradient by 


dc& = {g G X* \ a° (</)■, v) > (g,v ), for all v £ X}. 


In case a is differentiable this reduces to the standard gradient and in case cr is 
convex to the subgradient. 

In the literature on non-smooth optimization one can find several algorithms 
to minimize based on this generalized gradient. For some first numerical 

tests, we use a simple subgradient method: 

(j) k+1 =(t> k + h k A k A k £ d c a{^ k ) cr° = Id (28) 

where Id(x) = x is the identity transform. Note that this method does not use 
the full generalized gradient dccr but just one element of it in each step. This 
is typical for non-smooth/convex optimization methods because usually the full 
generalized gradient is not known. Since non-smooth optimization problems 
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often have kinks at the minimum itself, we cannot use standard techniques to 
control the step size and use a fixed rule 

hk = ak 0 , a > 0, 0 < /? < 1 (29) 

instead. This simple method converges for convex functions [2 |T4] and yields 
good results in the numerical experiments below. More sophisticated methods 
including convergence analysis for non-convex problems are available, see e.g. 

mm- 


5.2 Global minima? 

Because the objective function is non-convex in general, we must make 

sure that we do not end up in a suboptimal local minimum. In this section, we 
discuss some preliminary ideas to overcome this issue for a simple class of Id 
problems. To this end, let us consider piecewise constant functions 


Uo, 

for x < X\ (m) 


Ui, 

for Xi(fi) <x< x i+ i(^) 

(30) 

«n, 

for x n (n) < x 



with smooth parameter dependent jump locations Xi(fj,). We assume that the 
order Xq (/Li) < • • • < x n (fi) never changes and that the jump locations are well 
separated i.e. there is a constant L > 0 with \xi(/ji) — x-i- i(/x)| > L, i = 1,..., n. 


Transformed snapshot interpolation and training error Let us first 
state the transformed snapshot interpolation for these functions and the op¬ 
timization problem to find the inner transform. Because u{x 1 n) is piecewise 
constant in x, for transforms </>(/z, rf) that perfectly align the discontinuities the 
transformed snapshots v^x^rf) are constant in ry. Therefore, it is sufficient to 
confine ourselves to one single snapshot for the outer interpolation, say u(x, Mo) 
with node V n = {mo} so that we obtain 


u(x,n) = V^X^Lq) = u(((>(m,Mo),Mo)- 


(31) 


It follows that we have to compute inner transforms fJ. o) for all fj, G V and 

the single fixed node yto- To this end, we assume to know additional training 
snapshots «,(•, /ii),..., u(-, p. m ) for interpolation points /ii < ••• < /Lt m , with 
Mo ^ Mi f° r simplicity. Because we just use on snapshot for the outer interpo¬ 


lation the training error (24) reduces to 


<xr{4>)= sup sup ||u(-, Mi) — u (4>(ni, Mo); Mo) II Li (fi)- 


l<i<m 


l<i<m 


(32) 


Since all transforms o), 1 < 1 < m are uncorrelated, we can further 

simplify this and optimize for each transform </>(mo Mo) individually, which yields 

0(M;,Mo) = argmin||u(-,M l ) ~ u((p(-), Mo^L^n) (33) 
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for i = 1 Finally, with <j>(no, Ho){x) = x to ensure the interpolation 

condition ([5]), as in (16) we can define the full transform by an interpolation 


</Wm,Mo) = Ho) 


i —0 


where are the Lagrange polynomials for the nodes /z 0 ,..., /z n 


Counterexample: Local minima It remains to solve the m optimization 


problems (33). Already for this simple problem, optimization methods rely¬ 


ing on local search for updates can easily be fooled into a non-optimal local 
minimum. To this end, consider the example in Figure [4] defined by 

u(x,n) = x/i( M ) +X/ 2 (m)> j i(m) = [M,M + 1), = [m + 4,/x + 5), (34) 

where x is the characteristic function and the parameter shifts the entire func¬ 
tion. 


r 


Mo 

Mi 

M2 


Figure 4: Functions (34) for various parameters 


The snapshots /io and /z 2 in Figure @ are arranged such that the hrst 
interval of u(x,^ 2 ) intersects the second one of u(x,fio). Therefore the error 
er M2 (</>) is simply the area of the mismatch between the overlapping intervals 
plus the area of the two mismatched intervals. Note in particular that any 
small perturbation of </Hm2jMo) does not change the later error contribution. It 
follows that optimization schemes exclusively relying on local information are 
fooled into a local minimum that matches the (wrong) intersecting intervals. 

However, the situation changes, when the difference between the snapshot 
parameter /z q and the training parameter /it is small as e.g. for /zi in Figure 
|4j In this case there are no mismatched intervals and intuitively already simple 
subgradient based optimization schemes converge to the correct global minimum 
perfectly aligning the two functions. That this is in fact true is discussed in the 
following. 


Local convexity We confine ourselves to spatially monotone transforms Mo) 
in agreement with our assumption that the order of the jumps Xi(fi) does not 
change. Ideally, we search for transforms which exactly match the jump loca¬ 
tions 

<Mm>Mo )(xi(n)) = £?:(mo) ^ </>(m, Mo) -1 (ah(Mo)) = ^(m)- 
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Practically, we have to deal with perturbations so that the jumps only approx¬ 
imately match 

<KM>Mo) _ 1 (ab(Mo)) ~ 

If this matching error is sufficiently small compared to the minimal jump dis¬ 
tance L , such that only adjacent intervals of u(-,n) and o) overlap the 

training error simplifies to 

n 

= X! \ Ui ~ _ </KM,Mo) _1 (2h(Mo))|- (35) 


This error is convex in fi o) -1 which is not surprising because we assume 

that we are already close to a minimum. Whereas in principle the convexity 
allows us to compute optimal transforms this is not yet very practical since it 
requires very good initial values. However the identity transform id{x ) = x 
satisfied 

\xi(n) - fer 1 (£ i (/z 0 ))| = | Xi(n) - Xi(ti 0 )\ (36) 


which is sufficiently small to guarantee the error representation (351 for /r suffi¬ 
ciently close to fi o so that id can be used as an initial value in that case. 

Formally, in order to obtain a convex optimization problem, we augment the 
original training error minimization (33) with the following constraints 


a>(0) -> min 

0(/x,/Lto)^ 1 € C(H) strictly monotonically increasing (37) 

\xi(n o) - <KM,Mo) _1 (zi(Mo))| <B, i = 0,...,n 


for some constant B > 0. Both constraints are clearly convex in 0(/x,/xo) _1 - To 
show convexity of the objective function, note that by the triangle inequality 
we have 


I Xi(n) - </>(m,M o) 0 ))| < \xi{n) - Xi(n 0 )| + |a!i(^o) - </>(m.Mo) 1 (xi{n o))|- 

Thus, for [i sufficiently close to fj ,o and B sufficiently small the objective function 
reduces to (35) which is convex with respect to ^>(/x,Mo) -1 - According to (36) 
the identity is allowed by the constraints so that we know a suitable initial value 
for iterative optimization methods. Moreover for a transform </>(/z,^ 0 ) perfectly 
aligning the discontinuities, we have 


Xi{fi o) - <Km,Mo) (xi(fJ, o))| = \xi(n o) - Xi{n)\ 


which is also allowed by the constraints. It follows that the optimal error is 
cr M (^>) = 0 which is therefore a global minimum. 


Locality by transitivity In summary for ^ sufficiently close to /jq, we can 
reliably find a global minimum of the error a fl (<p) by solving a convex optimiza¬ 
tion problem, eventually with the identity id{x) = x as initial value. But what 
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about larger differences of /i and 7 x 0 as e.g. in our counter example with /i 2 in 
Figure |4]? To this end, recall that the main purpose of the transform 0 ( 74 , 77 ) is 
the alignment of jumps and kinks. Thus, if x(h) is the location of a jump for 
parameter 74 , we want the transforms to satisfy 

= x(rj). 

This condition guarantees that the transformed snapshot rj) = u{x{rj),n) 

has a jump at x(/j,) just as the target function u{x, 7 *). But this alignment con¬ 
dition is transitive in nature: For three consecutive parameters 7 / 0 , 7 x 1 and 7 x 2 
we have 

(0(/*i,Mo) 0 0(/*2,/*i))(z(/*2)) = x(no) 

so that 0 ( 741 , 7 * 0 ) o 0 ( 7 x 2,74 1 ) correctly aligns the jumps for parameters 7/0 and 
7 x 2 . Because this alignment property is a major requirement for the transform 
76 ( 7 x 2 , 7 * 0 ), we can define it that way 

0(7*2, Hi) := (^( 7 x 1 , 7 x 0 ) o 0(7*2, T*i)• ( 38 ) 


Let us apply this construction to our example transformed snapshot interpo¬ 


lation (311 where we have one snapshot at 7*0 and m training snapshots at 
7 x 1 < • • • < 7 i m . If we enforce transitivity (38 1 , we define 


0(7*2, 7 *o) := 0 (/*i,/*o) o • • • o 0 (t 4 , 7 *i-1) 

so that we are left with the calculation of the “local in 77 ” transforms 0 (Tq, Hi-i)- 
Because they are supposed to align jumps and kinks of and 

we can compute them by solving the optimization problem 


0 ( 7 * 2 , 7 * 2 - 1 ) = argmin ||i*(-, 74 ) - xx(<^(-), 7 x^ 1 ) 


l-Li(Q) 


(39) 


which is the same as the original problem (33) with 7 x 0 replaced by 7 x^— 1 . By 
choosing sufficiently many training snapshots, we can enforce 17 x 4 — 7 *x-i| to 


be sufficiently small such that the optimization problem (39) becomes convex. 


Therefore, we can reliably find global minimizers 0 ( 7 * 1 , 74 x_i) and in turn by 
(38) a transform that perfectly aligns the discontinuities for the parameters 7 x 0 
and Hi. This leads to a zero training error 07 -(0) which is therefore a global 
minimum. 

In summary, for the reconstruction of piecewise constant functions in Id 
from one snapshot, we can find 0 ( tx , 77) as the global minimum of the training 
error provided there are sufficiently many training snapshots. Of course the ar¬ 
gument relies on a couple of assumptions that are not true in more general cases. 
Notably, the functions u(x,n) might not be piecewise constant, we eventually 
want to use more snapshots for the outer interpolation and the parameter and 
spacial dimensions can be larger that one. Nonetheless, a transitivity property 
of the transforms is still realistic. As for the simple example of this section, this 
allows some locality in the parameter for the minimization of the training error. 
How to make use of this and to what extend this is helpful for more complicated 
scenarios is an open problem. 
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6 Numerical experiments 


In this section, we consider some first numerical tests of the transformed snap¬ 
shot interpolation ([7]). First, in Section 6.1 a Id example is presented, where 
the focus is on the approximation rate, while the inner transforms 77 ) are 
given explicitly. Then, in Section |6.2| the method is tested with a 2d Burgers 
Riemann problem where the solution is explicitly known. Finally, in Section [6.3| 
the method is applied to a shock bubble interaction which is a more challenging 
test case for the optimizer of the inner transform. 


6.1 Cut off Gaussian 

For a first numerical experiment, we consider the parametric function 


i(x, 77 ) := N 


0.4 + /r 


- 1 


N(x) := 


0.4 e 

0 


-7.0 a; 2 


— 1 < X < — \ 

else 


(40) 


which is a scaled and shifted Gaussian, cut off at a parameter dependent loca¬ 
tion, see Figure[5] This function is not chosen with a parametric PDE in mind, 
but has a parameter dependent jump and because it is known explicitly it is well 
suited for analyzing the performance of the transformed snapshot interpolation. 
For the snapshots we consider to alternatives: First we use the exact function 


(40) and second we interpolate it by piecewise linear functions on a uniform 


grid. The latter choice should simulate the outcome of PDE solvers which yield 
similar approximations of the parametric solution. Due to the simplicity of the 
example, we choose shifts for the inner transform: 

= x + s(n,r]). 


Recall that for the transformed snapshot interpolation (15), we only need to 
know s(/i, 77) for interpolation points 77 £ V m and 77 € V n . Thus, we can encode 
<f> by storing m x n floating point numbers. For all examples, we choose V n = V n . 
In addition, in this example we only consider the approximation properties of the 
transformed snapshot interpolation. In order not to interfere with the optimizer 
for actually finding the transform, we use an explicit formula for 5 ( 77 , 77 ) that 
exactly aligns the jumps and consider the optimizer in the numerical examples 
below. 

The numerical results are summarized in Figure[5] In case we use exact snap¬ 
shots, we see a more than polynomial convergence rate. Note that after aligning 
the snapshots, the transformed snapshots ^(*,77) are analytic in 77 so that this 
behaviour is in line with the error bounds of Proposition |2.1[ However, for the 
linearly interpolated snapshots the situation is different. The error first decays 
and then saturates at a level dependent on the spacial grid resolution. These 
levels correspond to the maximal error of the snapshots themselves as shown 
in Table [TJ This makes sense because the transformed snapshot interpolation 
error can hardly be better than the error of the snapshots it relies on. 
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For a comparison, Figure [5] also contains the error of a simple polynomial in¬ 
terpolation without transform. We see the typical staircasing behaviour and an 
error that is orders of magnitudes worse than the one with previous transform. 


•1CT 2 




Figure 5: Left: Errors of transformed snapshot interpolations of (40). Black 
with dot marks: Snapshots are piecewise linear interpolations with uniform grid 
sizes h = 0.1,0.01,0.001,0.0001. Red with + marks: Snapshots are exact func¬ 
tions without any approximation. Blue with x marks: Interpolation without 
transform. Right: Truth solution (black) and interpolation from nine snapshots 
with (red, dashed) and without (blue dotted) transform. 



transformed snapshot interpolation 

interpolation 

n 

0.1 

0.01 

0.001 

0.0001 

exact 

exact 

2 

9.4- 10" 4 

3.83 • 10" 4 

3.81 ■ 10" 4 

3.81 ■ 10" 4 

3.61 • KT 4 

1.02 • 10" 2 

3 

3.6 • 10 -3 

3.19- 10" 4 

4.18 • 10~ 5 

4.17- 10" 5 

3.3 • 10“ 5 

7.47- 10" 3 

5 

5.92 • 10” 3 

5.13- 10" 4 

2.96 • 10~ 5 

2.31 • 10" 6 

1.49 • KT 7 

2.2 • KT 3 

9 

3.43 • 10” 3 

3.74- 10" 4 

3.77-10” 5 

3•10" 6 

1.07- 10~ n 

2.51 • 10" 3 


maximal L i(fl) error of the snapshots 



3.58 • KT 3 

3.53 • 10 -4 

3.48 • 10” b 

3.66 ■ KT iU 




Table 1: Errors for example (40) for different uniform spacial grids and num¬ 
ber of snapshots (n). The last line contains the maximal Li(fl) error of the 
respective snapshots. 


6.2 2d Burgers’ equation 

For a second example, we consider the two dimensional Burgers’ equation 

d t u + V • f J = 0 
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with v = ( 1 , 1 ) T and the initial condition 




on the unit cube O = [ 0 , l] 2 . 
problem is 


- 0.2 

if x 

< 

0.5 

and 

y 

> 

0.5 

- 1.0 

if x 

> 

0.5 

and 

y 

> 

0.5 

0.5 

if x 

< 

0.5 

and 

y 

< 

0.5 

0.8 

if x 

> 

0.5 

and 

y 

< 

0.5 

According 

to 

m 

IB] the 

s exact 


= < 


-0.2 

0.5 

- 1.0 

0.5 

- 1.0 

0.5 

- 1.0 

2x -1 
24 

- 1.0 

0.8 


if x < | — f 


if \ 


y < X < 


and 


— 4 and 


V > 2 


if I 


i 1 < a; < 4 


and 


3t 
20 : 

otherwise, 

y>-f + 

otherwise, 

^ _L A 
y ^ 6 ^ 12 

otherwise, 


15t 
28 5 


5 1 
24’ 


if I 


4+4<x<}+y and 


V > x 184 ^ 2 ) 

otherwise, 


if 5 + f < * 


and 


y > £ - 


_t_ 

10 ’ 
otherwise, 


(41) 

For a simple test of the transformed snapshot interpolation, we consider the time 
t as the parameter of interest so that the snapshots are solutions at various time 
instances used for the reconstruction of the solution at intermediate times. Note 
that with this choice of the parameter the solution has exactly the features in 
question: it has parameter dependent jumps and kinks along non-trivial curves. 
In addition the exact solution is known which is helpful for an exact assessment 
of the numerical errors. In order to simulate a numerical solution of the Burger’s 
equation, we sample the snapshots on a 100 x 100 grid and use a piecewise linear 
reconstruction from these samples. Also the integrals for evaluating the errors 
during the optimization of the inner transform <j> rely on this grid. 

Figure [ 6 ] shows the results for a reconstruction at time 0.45 from two snap¬ 
shots at times 0.3 and 0.5 with an additional training snapshot at time 0.4 to 
define the training error For a first test, the inner transforms <p(n,ri) 

for fi,r] £ {0.3, 0.5} are simply polynomials mapping R 2 —> R 2 . In general 
this choice does not guarantee that fl is mapped to itself, however it is easy to 
enforce that the edges of the rectangular domain are mapped to itself so that 
small perturbations of the identity are diffeomorphisms. In order to be able to 
align both the kink and the jump in the lower right corner, for the x-component 
we choose a (multivariate) polynomial of degree 3x2 and for the y- com P onen l 
a polynomial of degree 2x2. For the optimization of the training error with 


respect (j) we use a subgradient method (28 1 , (29) with 500 steps of the rather 
conservative fixed step size 

; _ 10" 3 
k ~ (* + l)0-i ■ 

Compared to a classical polynomial interpolation, the additional transform al¬ 
most completely removes the artificial staircasing behaviour. Also the kinks 
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around the “ramp” in the upper left corner of the figures are much better re¬ 
solved. The L i errors computed by an adaptive quadrature instead of the grid 
of the snapshots are as follows. 

L|-error interpolation 0.0355373675439 

L|-error transformed snapshot interpolation 0.00739513000396 
maximal snapshot Li-error 0.0051148730424 


In conclusion the additional inner transform reduces the error almost by a fac¬ 
tor of five compared to a plain polynomial interpolation. Note that the error 
of the transformed snapshot interpolation is almost down to the maximal error 
of the snapshots themselves. As we have verified in Figure [5] for the Id exam¬ 
ple of Section 6.1 we expect the error to saturate somewhere around this level 
so that more snapshots or degrees of freedom for the inner transform are not 
expected to yield major improvements. This is a serious bottleneck for com¬ 
puting convergence rates: Due to the jump discontinuities the maximal error of 
the snapshots converges with a low rate. The resulting high number of spacial 
degrees of freedom renders the computation of convergences rates challenging. 



Figure 6: Reconstruction of the solution (411 for Burgers equation at the time 
t = 0.45. Left: polynomial interpolation, middle: transformed snapshot inter¬ 
polation, right: exact solution. 


6.3 Shock bubble interaction 

For a last more complicated example, we consider a compressible Euler simula¬ 
tion of a shock-bubble interaction m ■ Because the code for the above examples 
relies on piecewise linear interpolation on a uniform grid to represent the snap¬ 
shots it is straight forward to read them from pictures. For the shock bubble 
interaction experiments they are frames from a video showing the time evolu¬ 
tion of the density provided by [20] [T9] . Figure [7] shows the snapshots and the 
reconstruction at a new time by linear interpolation and transformed snapshot 
interpolation. Comparable to Section |6.2[ we simply choose third order poly¬ 
nomials for the inner transform <f> mapping the edges of the domain to itself. 
We see that the linear interpolation result basically shows the two bubbles from 
the original snapshots, whereas the true solution of course just has one bubble. 
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Using the additional transform </>, the second picture finds the correct location 
of the shock and the bubble. Thus despite lots of more complicated fine struc¬ 
ture in the pictures, the optimizer reliably finds the correct transform. Having 
a closer look, the reconstructed bubble appears to be a little blurred. However, 
we only use third order polynomials which eventually is insufficient for a perfect 
alignment of the shapes. 



Figure 7: Left column: Snapshots for time indices 3.5, 3.75 and 4.0 where the 
middle one is only used for the optimization of the transform </>. Right column: 
reconstruction at time index 3.7 by linear interpolation (top) and transformed 
snapshot interpolation (bottom). 


Appendix A Linear Width 

As an example for the limitations of polyadic decomposition based methods, let 
us consider their best possible performance for the following simple parametric 
transport problem: 

:= Ut + gu x =0 for 0 < t < 1, x € K. 

u{x, 0) = g(x) : = | x X x > o for t = 0, 

with parameter g, G V = [/x m i n , £i max ] C R. Its solution is given by 

u(x,t) = g(x — gt). (42) 

The typical benchmark for the performance of reduced basis methods is the 
Kolmogorov n-width 

d n {F) = inf sup inf ||u - (j>\\ Ll 

dim Y—n U (zj? <f)£Y 

of the solution manifold 

T := M-,/r )\g GV} = {g(x - gt)\g G V}. (43) 
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However, with X n := span{^>j : i = 1,... ,n} based on the ipi of the polyadic 
decomposition ([I]), we conclude that 


d n (F) < sup inf 

uer<l>ex n 


IK-,m) - 


\ Ll < sup 


*(■>/*) 


Therefore, the errors of polyadic decomposition based methods, including but 
not restricted to reduced basis methods, are lower bounded by the Kolmogorov 
ro-width if the error is measured in the sup-nornr with respect to the parameter 
variable. Of course this measure for the error is not appropriate for all applica¬ 
tions, but we use it here to exemplify the limitations of the standard polyadic 
decomposition based methods. 

For our simple model problem the Kolmogorov n-width is bounded according 
to the following Proposition. The proof is similar to ,7], see also (TS]. Note 
that order one is already achieved by a simple nonadaptive piecewise constant 
approximation as e.g. in ©>■ 


Proposition A.l. The Kolmogorov width of the solution manifold T defined 

-l 


in (43) satisfies 


d n {F) 


n 


i.e. d n {F) is equivalent to n 1 up to a constant. 


Proof. We show the lower bound by comparing the Kolmogorov n-width of the 
solution manifold J- to the known width of a ball in the Li-norm. For the 
construction of this ball, we choose a uniform distribution of k + 1 snapshots 
with parameters 


— /^min 


+ k { ^ 


Mmin)i i 0, . . . , /t, 


where k will be chosen below. Note that these specific snapshots are only used 
for this proof and are not necessarily used in actual approximation methods like 
e.g. reduced basis methods. We use the snapshots to define the functions 


6 - ’UfJ’i ^Ui-1 5 ^ - 1, . . . , fc, 

which will be the corners of the Li-ball. From 


inf ||6 - : 

y£Y 


L i ^ 


< inf I 

y£Y 




\ Ll + inf ||w 

y&Y 


V'i-i 


y\\ Ll < 2 sup inf ||u — y\\L 1 

u^TV^Y 


for any linear space Y of dimension at most n follows that 

d„({£o, • ■ • ,£k}) < ^d n {T). 

We complete the set {£o,..., £&} to a full ball without increasing the Kolmogorov 
width. To this end assume that A*, i = 1,..., k satisfy |Ai| < 1 and yi , 


26 







i = 1,..., k are the minimizers of vnty^y ||6 — y|ki- Then we have 


inf 

y&Y 




< 


k k 

E A *& " E A '^' 

2=1 2=1 


<^1^1116-^11^ <d„({6,---,6J) < 2 d n (B) 


for any space Y of dimension smaller than n realizing the Kolmogorov width of 
T. It follows that for the set 

! fc k 

E A *& : E i A ii < 1 

2=1 2=1 


we have 

dniF) < 2 d n (T). 


(44) 


Next, we show that T' is in fact an Li-ball. To this end, note that according to 
the exact solution (421 the functions 6 take the value one on disjoint triangles 
and zero else. The area of the triangles is h/2 with h = (/z max —/i m i n )/k, because 
the time t is bounded between 0 < t < 1. Thus, we have 


1161 k, 


h 

2 ' 


It follows that 


E A *k 


2=1 


= El«k = ~EM- 


2=1 


2=1 


so that 


T' = span{6, • • •, 6J n B h L [ 2 


where B E is the Li-ball with radius h/2. Choosing k = 2 n, we obtain the 
Kolmogorov width 

dn(T') = 


see e.g. m- Using ( |44[ ) and h ~ l/n completes the proof of the lower bound. 

In order to prove the upper bound, note that — u is one on a triangular 
domain and zero else where u tJi are the snapshots used in the proof of the lower 
bound. Calculating the area of the triangle as before, this yields 


u u - u u 


h 

Li < 77 


for the parameter fii closest to /r. Thus, a piecewise constant approximation by 
..., Ufj, k yields the upper bound of the proposition. ;P 
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